---
title: "R Notebook"
output: html_notebook
---
#Load packages
```{r}
library(ggthemes)
library(lme4)
library(optimx)
library(tidyverse)
library(emmeans)
```

#Load and Wrangle Data

##General Wrangling
```{r}
Data <- read.csv("../Samson_data_final.csv")

#Prepare general DS
Data <- Data %>% 
  mutate_at( #Ensure proper data type for factors
    .vars = vars(Sex, Subject, PerspectiveSelf), factor) %>% 
  add_column( #Create a unique studyid
    PaperStudy = paste(Data$Paper, Data$StudyNumber, sep="_")) %>% 
  mutate( #re-code consistency levels as numeric, for use in analyses
    Directional.Consistency = recode(Directional.Consistency, 'I' = -1, 'N' = 0, 'C' = 1),
    Perspective.Consistency = recode(Perspective.Consistency, 'I' = -1, 'N' = 0, 'C' = 1),
  )

#Prepare a DS which contains only RT data, and apply exclusions for analyses
DataRT <- Data %>% 
  filter(
    RTorError=="RT",
    PerspectiveSelf==1,
    Neurotypical=="Y")

#Prepare a DS which contains only ER data, and apply exclusions for analyses
DataError <- Data %>% 
  filter(
    RTorError=="Error",
    Value < 1,
    PerspectiveSelf==1,
    Neurotypical=="Y")
```


##Filtering for "Filtered Analyses"
```{r}

#for each case of persp consistency, are there multiple instances of dir consistency? record these cases
varied_dir_RT <- DataRT %>% 
  group_by(Paper, StudyNumber, Perspective.Consistency) %>%
  summarise(
    dir_values = n_distinct(Directional.Consistency)) %>% 
  mutate(filter_condition = paste(Paper, StudyNumber, sep="_")) %>% 
  filter(dir_values > 1)

#for each case of dir consistency, are there multiple instances of persp consistency? record these cases
varied_persp_RT <- DataRT %>% 
  group_by(Paper, StudyNumber, Directional.Consistency) %>%
  summarise(
    persp_values = n_distinct(Perspective.Consistency)) %>% 
  mutate(filter_condition = paste(Paper, StudyNumber, sep="_")) %>% 
  filter(persp_values > 1)

#for each case of persp consistency, are there multiple instances of dir consistency? record these cases
varied_dir_Error <- DataError %>% 
  group_by(Paper, StudyNumber, Perspective.Consistency) %>%
  summarise(
    dir_values = n_distinct(Directional.Consistency)) %>% 
  mutate(filter_condition = paste(Paper, StudyNumber, sep="_")) %>% 
  filter(dir_values > 1)

#for each case of dir consistency, are there multiple instances of persp consistency? record these cases
varied_persp_Error <- DataError %>% 
  group_by(Paper, StudyNumber, Directional.Consistency) %>%
  summarise(
    persp_values = n_distinct(Perspective.Consistency)) %>% 
  mutate(filter_condition = paste(Paper, StudyNumber, sep="_")) %>% 
  filter(persp_values > 1)


#Create data sets for analyses without no direction/perspective
simpleRT <- DataRT %>% filter(Directional.Consistency!=0) %>% filter(Perspective.Consistency!=0) %>%
            mutate(Directional.Consistency=factor(Directional.Consistency),
                   Perspective.Consistency=factor(Perspective.Consistency))

simpleError <- DataError %>% filter(Directional.Consistency!=0) %>% filter(Perspective.Consistency!=0) %>%
            mutate(Directional.Consistency=factor(Directional.Consistency),
                   Perspective.Consistency=factor(Perspective.Consistency))

#Create data sets for filtered analyses (simply remove all cases where PC and DC covary)
varied_DataError <- DataError %>% 
  mutate(filter_condition = paste(Paper, StudyNumber, sep="_")) %>% 
  filter(filter_condition %in% varied_dir_Error$filter_condition | filter_condition %in% varied_persp_Error$filter_condition)

varied_DataRT <- DataRT %>% 
  mutate(filter_condition = paste(Paper, StudyNumber, sep="_")) %>% 
  filter(filter_condition %in% varied_dir_RT$filter_condition | filter_condition %in% varied_persp_RT$filter_condition)
```


#Analyses

##Summary Info
```{r}
DataSummary <- Data %>% 
  mutate(Experiment.Unique = paste(Paper, StudyNumber, '_'))

#number of papers
n_distinct(DataSummary$Paper)

#number of experiments
n_distinct(DataSummary$Experiment.Unique)

#number of subjects
n_distinct(DataSummary$Subject)

#papers with Error data
n_distinct(DataSummary[DataSummary$RTorError=='Error', "Paper"])

#papers with RT data
n_distinct(DataSummary[DataSummary$RTorError=='RT', "Paper"])
unique(DataSummary[DataSummary$RTorError=='RT', "Paper"])
```

##Unfiltered Analyses
  
###Response Time

####lm models

```{r}
#lm with Perspective and Direction
lm.1 <- lmer(Value ~ Perspective.Consistency +Directional.Consistency+ (Perspective.Consistency + Directional.Consistency |PaperStudy) + (1|Subject), data=DataRT,control=lmerControl(optimizer="bobyqa"))

#lm with only Perspective
lm.2 <- lmer(Value ~ Perspective.Consistency + (Perspective.Consistency + Directional.Consistency |PaperStudy) + (1|Subject), data=DataRT,control=lmerControl(optimizer="bobyqa"))

#lm with only Direction
lm.3 <- lmer(Value ~ Directional.Consistency + (Perspective.Consistency + Directional.Consistency |PaperStudy) + (1|Subject), data=DataRT,control=lmerControl(optimizer="bobyqa"))
```

####Anovas
```{r}
anova(lm.1, lm.2) #effect of directional consistency above and beyond perspective 
anova(lm.1, lm.3) #effect of perspective consistency above and beyond direction
```


###Error Rate
####lm models
```{r}
#lm with Perspective and Direction
lm.4 <- lmer(Value ~ Perspective.Consistency + Directional.Consistency+ (Perspective.Consistency + Directional.Consistency |PaperStudy) + (1|Subject), data=DataError,control=lmerControl(optimizer="bobyqa"))

#lm with only perspective
lm.5 <- lmer(Value ~ Perspective.Consistency + (Perspective.Consistency + Directional.Consistency |PaperStudy) + (1|Subject), data=DataError, control=lmerControl(optimizer="bobyqa"))

#lm with only direction
lm.6 <- lmer(Value ~ Directional.Consistency + (Perspective.Consistency + Directional.Consistency |PaperStudy) + (1|Subject), data=DataError, control=lmerControl(optimizer="bobyqa"))
```

####Anovas
```{r}
anova(lm.4, lm.5) #effect of directional consistency above and beyond perspective 
anova(lm.4, lm.6) #effect of perspective consistency above and beyond direction
```

## Simple Analyses

###Response Time

####lm models

```{r}
#lm with Perspective and Direction
lm.7 <- lmer(Value ~ Perspective.Consistency +Directional.Consistency+ (Perspective.Consistency + Directional.Consistency |PaperStudy) + (1|Subject), data=simpleRT,control=lmerControl(optimizer="bobyqa"))

#lm with only Perspective
lm.8 <- lmer(Value ~ Perspective.Consistency + (Perspective.Consistency + Directional.Consistency |PaperStudy) + (1|Subject), data=simpleRT,control=lmerControl(optimizer="bobyqa"))

#lm with only Direction
lm.9 <- lmer(Value ~ Directional.Consistency + (Perspective.Consistency + Directional.Consistency |PaperStudy) + (1|Subject), data=simpleRT,control=lmerControl(optimizer="bobyqa"))
```

####Anovas
```{r}
anova(lm.7, lm.8) #effect of directional consistency above and beyond perspective 
anova(lm.7, lm.9) #effect of perspective consistency above and beyond direction
```

###Error Rate
####lm models
```{r}
#lm with Perspective and Direction
lm.10 <- lmer(Value ~ Perspective.Consistency + Directional.Consistency + (Perspective.Consistency + Directional.Consistency |PaperStudy) + (1|Subject), data=simpleError,control=lmerControl(optimizer="bobyqa"))

#lm with only perspective
lm.11 <- lmer(Value ~ Perspective.Consistency + (Perspective.Consistency + Directional.Consistency |PaperStudy) + (1|Subject), data=simpleError, control=lmerControl(optimizer="bobyqa"))

#lm with only direction
lm.12 <- lmer(Value ~ Directional.Consistency + (Perspective.Consistency + Directional.Consistency |PaperStudy) + (1|Subject), data=simpleError, control=lmerControl(optimizer="bobyqa"))
```

####Anovas
```{r}
anova(lm.10, lm.11) #effect of directional consistency above and beyond perspective 
anova(lm.10, lm.12) #effect of perspective consistency above and beyond direction
```

##Filtered Analyses

###Response Time
```{r}
#lm with Perspective and Directional
lm.13 <- lmer(Value ~ Perspective.Consistency + Directional.Consistency + (Perspective.Consistency + Directional.Consistency |PaperStudy) + (1|Subject), data=varied_DataRT, control=lmerControl(optimizer="bobyqa"))

#lm with only Perspective
lm.14 <- lmer(Value ~ Perspective.Consistency + (Perspective.Consistency + Directional.Consistency |PaperStudy) + (1|Subject), data=varied_DataRT, control=lmerControl(optimizer="bobyqa"))

#lm with only Directional
lm.15 <- lmer(Value ~ Directional.Consistency + (Perspective.Consistency + Directional.Consistency |PaperStudy) + (1|Subject), data=varied_DataRT, control=lmerControl(optimizer="bobyqa"))

```

####Anovas
```{r}
anova(lm.13, lm.14) #effect of directional consistency above and beyond perspective 
anova(lm.13, lm.15) #effect of perspective consistency above and beyond direction
```

###Error Rate
```{r}
#lm with Perspective and Directional
lm.16 <- lmer(Value ~ Perspective.Consistency + Directional.Consistency+ (Perspective.Consistency + Directional.Consistency |PaperStudy) + (1|Subject), data=varied_DataError, control=lmerControl(optimizer="bobyqa"))

#lm with only perspective
lm.17 <- lmer(Value ~ Perspective.Consistency + (Perspective.Consistency + Directional.Consistency |PaperStudy) + (1|Subject), data=varied_DataError, control=lmerControl(optimizer="bobyqa"))

#lm with only directional
lm.18 <- lmer(Value ~ Directional.Consistency + (Perspective.Consistency + Directional.Consistency |PaperStudy) + (1|Subject), data=varied_DataError, control=lmerControl(optimizer="bobyqa"))
```

####Anovas
```{r}
anova(lm.16, lm.17) #effect of directional consistency above and beyond perspective 
anova(lm.16, lm.18) #effect of perspective consistency above and beyond directionLM
```

##Automatic Theory of Mind Analysis

###RT
```{r}
## RT effect
lm.19 <- lmer(Value ~ Perspective.Consistency * Explicitly.Tracking.Other +  (Perspective.Consistency + Explicitly.Tracking.Other |PaperStudy) + (1|Subject), data=DataRT, control=lmerControl(optimizer="bobyqa"))

## Interaction effect
lm.20 <- lmer(Value ~ Perspective.Consistency + Explicitly.Tracking.Other +  (Perspective.Consistency + Explicitly.Tracking.Other |PaperStudy) + (1|Subject), data=DataRT, control=lmerControl(optimizer="bobyqa"))

anova(lm.19,lm.20)

##Factor version for pairwise comparison
lm.19f <- lmer(Value ~ factor(Perspective.Consistency) * Explicitly.Tracking.Other +  (Perspective.Consistency + Explicitly.Tracking.Other |PaperStudy) + (1|Subject), data=DataRT, control=lmerControl(optimizer="bobyqa"))

##NB: We simplified the random effects structure by removing the interaction slope from the PaperStudy effect to help with convergence issues. However, the story is quite similar the convergence warnings if you include the random slope for the interaction.

# pairwise comparisons to decompose interaction:
emmip(lm.19f, Explicitly.Tracking.Other ~  Perspective.Consistency)
emmeans(lm.19f, pairwise ~ Perspective.Consistency | Explicitly.Tracking.Other, pbkrtest.limit = 12000)

```


### Error
```{r}
lm.21 <- lmer(Value ~ Perspective.Consistency * Explicitly.Tracking.Other +  (Perspective.Consistency + Explicitly.Tracking.Other |PaperStudy) + (1|Subject), data=DataError, control=lmerControl(optimizer="bobyqa"))

# Interaction effect
lm.22 <- lmer(Value ~ Perspective.Consistency + Explicitly.Tracking.Other +  (Perspective.Consistency + Explicitly.Tracking.Other |PaperStudy) + (1|Subject), data=DataError, control=lmerControl(optimizer="bobyqa"))

anova(lm.21,lm.22)

## Factor version for pairwise comparison
lm.21f <- lmer(Value ~ factor(Perspective.Consistency) * Explicitly.Tracking.Other +  (Perspective.Consistency + Explicitly.Tracking.Other |PaperStudy) + (1|Subject), data=DataError, control=lmerControl(optimizer="bobyqa"))

# pairwise comparisons to decompose interaction:
emmip(lm.21f, Explicitly.Tracking.Other ~  Perspective.Consistency)
emmeans(lm.21f, pairwise ~ Perspective.Consistency | Explicitly.Tracking.Other, pbkrtest.limit = 12000)
```


#Figures
##Simple Data Visualization -- Visualize all data, excluding 'No perspective/direction'

###RT 
```{r}
## lm.7 is the full model for RT in the simplified analyses
## we use these emmeans estimates rather than the simple means to plot the lines which show the model estimated effects over the data+violins
#meansRT_estDir <- emmeans(lm.7, ~ Directional.Consistency,lmerTest.limit=10000,pbkrtest.limit = 10000)
#meansRT_estPer <- emmeans(lm.7, ~ Perspective.Consistency,lmerTest.limit=10000,pbkrtest.limit = 10000)


## order here is Consistent-Directional, Consistent-Perspective, Inconsistent-Directional, Inconsistent-Perspective
## based upon the above emmeans()
mod.estimates_rt=c(672,682,701,691)


#Prettier names for graph
DataRT_fig <- DataRT %>% 
  pivot_longer(cols = c(Perspective.Consistency, Directional.Consistency), names_to = 'Classification', values_to = 'Consistency') %>% 
  filter(Consistency != 0) %>% 
  mutate(
    Classification = case_when(
      Classification == "Perspective.Consistency" ~ "Perspective Consistency",
      Classification == "Directional.Consistency" ~ "Directional Consistency"),
    Consistency = case_when(
      Consistency == 1 ~ 'Consistent',
      Consistency == -1 ~ 'Inconsistent')) 

#calculate means 
RT_means <- DataRT_fig %>%
  group_by(Consistency, Classification) %>%
  dplyr::summarize(mean = mean(Value)) %>%
  ungroup()%>%
  mutate(mod.estimates=mod.estimates_rt)



#plot
plot1 <- ggplot(
  data = DataRT_fig,
  aes(x = factor(Consistency),
      y = Value,
      color = factor(Consistency))) + 
  geom_point(
    position = position_jitterdodge(dodge.width = 0.9, jitter.width = 0.9),
    size = 0.01,
    alpha = 0.5) +
  geom_violin(
    aes(alpha = 0),
    position = position_dodge(width=0.9)) +
  geom_boxplot(
    aes(alpha = 0),
    outlier.shape = NA,
    width = 0.15,
    position = position_dodge(width=0.9),
    fatten = NULL) +
  geom_line(
    data = RT_means,aes(x = Consistency, 
      y = mod.estimates, 
      group=Classification),
    color="grey30"
    ) +
  geom_point(
    data = RT_means, 
    mapping = aes(
      x = Consistency, 
      y = mod.estimates, 
      color = factor(Consistency)),
    position = position_dodge(width=0.9)
  ) +
  facet_grid(~Classification) +
  ylim(250, 1200) +
  theme(legend.position = "none") +
  xlab(element_blank()) +
  ylab("Response Time (ms)") +
  theme_few() +
  scale_colour_ptol() +
  theme(legend.position = "none")

plot1

#ggsave(filename = "../Figures/RT_simple_analysis.pdf", width = 5, height = 5)
```

###Error
```{r}
## lm.10 is the full model for Error in the simplified analyses
## we use these emmeans estimates rather than the simple means to plot the lines which show the model estimated effects over the data+violins
#meansER_estDir <- emmeans(lm.10, ~ Directional.Consistency,lmerTest.limit=10000,pbkrtest.limit = 10000)
#meansER_estPer <- emmeans(lm.10, ~ Perspective.Consistency,lmerTest.limit=10000,pbkrtest.limit = 10000)


## order here is Consistent-Directional, Consistent-Perspective, Inconsistent-Directional, Inconsistent-Perspective
## based on above emmeans()
mod.estimates_er=c(0.0545,0.0603,0.0816,0.0757)


#Prettier names for graph
DataError_fig <- DataError %>% 
  pivot_longer(cols = c(Perspective.Consistency, Directional.Consistency), names_to = 'Classification', values_to = 'Consistency') %>% 
  filter(Consistency != 0) %>% 
  mutate(
    Classification = case_when(
      Classification == "Perspective.Consistency" ~ "Perspective Consistency",
      Classification == "Directional.Consistency" ~ "Directional Consistency"),
    Consistency = case_when(
      Consistency == 1 ~ 'Consistent',
      Consistency == -1 ~ 'Inconsistent')) 

#calculate means
Error_means <- DataError_fig %>% 
  group_by(Consistency, Classification) %>% 
  dplyr::summarize(mean = mean(Value)) %>%
  ungroup() %>%
  mutate(mod.estimates=mod.estimates_er)


#plot
plot2 <- ggplot(
  data = DataError_fig,
  aes(x = factor(Consistency),
      y = Value,
      color = factor(Consistency))) + 
  geom_point(
    position = position_jitterdodge(dodge.width = 0.9, jitter.width = 0.6),
    size = 0.01,
    alpha = 0.5) +
  geom_violin(
    aes(alpha = 0),
    position = position_dodge(width=0.9)) +
  geom_boxplot(
    aes(alpha = 0),
    outlier.shape = NA,
    width = 0.15,
    position = position_dodge(width=0.9),
    fatten = NULL) +
  geom_line(
    data = Error_means,aes(x = Consistency, 
      y = mod.estimates_er, 
      group=Classification),
    color="grey30"
    ) +
  geom_point(
    data = Error_means, 
    mapping = aes(
      x = Consistency, 
      y = mod.estimates_er, 
      color = factor(Consistency)),
    position = position_dodge(width=0.9)
  ) +
  facet_grid(~Classification) +
  ylim(0, 0.5) +
  theme_few() +
  scale_colour_ptol() +
  theme(legend.position = "none") +
  ylab("Error Rate") +
  xlab(element_blank())


plot2

#ggsave(filename = "../Figures/Error_simple_analysis.pdf", width = 5, height = 5)
```

## Visualize the effect of 'explicitly tracking other'

### RT
```{r}

#Rename and reorder for pretty graph
autoToM_RT <- DataRT %>%
  mutate(
    Perspective.Consistency = factor(case_when(
        Perspective.Consistency == 1 ~ 'Consistent',
        Perspective.Consistency == 0 ~ 'No Perspective',
        Perspective.Consistency == -1 ~ 'Inconsistent')),
    Explicitly.Tracking.Other = factor(case_when(
      Explicitly.Tracking.Other == 'N' ~ "No Explicit Tracking",
      Explicitly.Tracking.Other == 'Y' ~ "Explicit Tracking")),
    Explicitly.Tracking.Other = fct_relevel(Explicitly.Tracking.Other, "No Explicit Tracking", "Explicit Tracking")
  ) %>% 
  filter(Perspective.Consistency != 'No Perspective')

#gather summary stats for plot
autoTom_sumRT <- autoToM_RT  %>% 
  group_by(Explicitly.Tracking.Other, Perspective.Consistency) %>%
  summarize(
    n=length(Value),
    mean = mean(Value),
    sd = sd(Value),
    se = sd/sqrt(n))

#create plot
plot3 <- autoToM_RT %>% 
  ggplot(
    aes(
      x=factor(Perspective.Consistency),
      y=Value,
      color=Perspective.Consistency)) +
  geom_violin() +
  geom_point(
    size = 0.2,
    position = position_jitterdodge(
        jitter.width = .9,
        dodge.width = .9),
        alpha=.18) +
  geom_boxplot(
    aes(alpha = 0),
    outlier.shape = NA,
    width = 0.15,
    position = position_dodge(width=0.9),
    fatten = NULL) +  
  geom_point(
    data=autoTom_sumRT, 
    aes(
      y=mean,
      color = factor(Perspective.Consistency),
      x=factor(Perspective.Consistency)),
    size=.8,
    position=position_dodge(width=.9)) +
  geom_line(
    data=autoTom_sumRT,
    color="grey30",
    aes(
      y=mean,
      group=Explicitly.Tracking.Other),
    size=0.5,
    position = position_dodge(width = .9)) +
  facet_grid(~Explicitly.Tracking.Other) +
  ylab("Response Time") +
  xlab(element_blank()) +
  ylim(250, 1500) +
  theme_few() +
  scale_colour_ptol() +
  theme(legend.position = "none") 
  
plot3

#ggsave(filename = "../Figures/RT_tracking_other.pdf", width = 5, height = 5)
```

### Error
```{r}

#Rename and reorder for pretty graph
autoToM_Error <- DataError %>%
  mutate(
    Perspective.Consistency = factor(case_when(
        Perspective.Consistency == 1 ~ 'Consistent',
        Perspective.Consistency == 0 ~ 'No Perspective',
        Perspective.Consistency == -1 ~ 'Inconsistent')),
    Explicitly.Tracking.Other = factor(case_when(
      Explicitly.Tracking.Other == 'N' ~ "No Explicit Tracking",
      Explicitly.Tracking.Other == 'Y' ~ "Explicit Tracking")),
    Explicitly.Tracking.Other = fct_relevel(Explicitly.Tracking.Other, "No Explicit Tracking", "Explicit Tracking")
  )%>% 
  filter(Perspective.Consistency != 'No Perspective')

#gather summary stats for plot
autoTom_sum_Error <- autoToM_Error  %>% 
  group_by(Explicitly.Tracking.Other, Perspective.Consistency) %>%
  summarize(
    n=length(Value),
    mean = mean(Value),
    sd = sd(Value),
    se = sd/sqrt(n))

#create plot
plot4 <- autoToM_Error %>% 
  ggplot(
    aes(
      x=factor(Perspective.Consistency),
      y=Value,
      color=Perspective.Consistency)) +
  geom_violin() +
  geom_point(
    size = 0.2,
    position = position_jitterdodge(
        jitter.width = .9,
        dodge.width = .9),
        alpha=.18) +
  geom_boxplot(
    aes(alpha = 0),
    outlier.shape = NA,
    width = 0.15,
    position = position_dodge(width=0.9),
    fatten = NULL) +  
  geom_point(
    data=autoTom_sum_Error, 
    aes(
      y=mean,
      color = factor(Perspective.Consistency),
      x=factor(Perspective.Consistency)),
    size=.8,
    position=position_dodge(width=.9)) +
  geom_line(
    data=autoTom_sum_Error,
    color="grey30",
    aes(
      y=mean,
      group=Explicitly.Tracking.Other),
    size=0.5,
    position = position_dodge(width = .9)) +
  ylab("Error Rate") +
  xlab(element_blank()) +
  ylim(0, .5) +
  labs(color = "Explicitly Tracking Other") +
  facet_grid(~Explicitly.Tracking.Other) +
  theme_few() +
  scale_colour_ptol() +
  theme(legend.position = 'none')
  
plot4

#ggsave(filename = "../Figures/Error_tracking_other.pdf", width = 5, height = 5)
```